Introducción: Análisis previo de los datos

En este apartado nos centraremos en entender la composición de nuestros datos, así como la tendencia, estacionariedad y estacionalidad de los mismos, con el fin de ajustar posteriormente diversos modelos predictivos para finalmente quedarnos con el más efectivo a la hora de predecir. Por el camino realizaremos algunos test estadísticos y gráficos que nos ayudarán a determinar cuál debe ser el tratamiento de la serie, prevaleciendo siempre el criterio de minimizar el error de predicción.

## Warning: package 'ggplot2' was built under R version 3.5.2

Como vemos en el gráfico, la tendencia es claramente descendente, mostrando que en los últimos años el uso del cine como entretenimiento ha decaído considerablemente. Una limitación de esta serie temporal es que llega sólo hasta el año 2012 (este año está completo), justo el año en el que entraron empresas como Netflix al mercado. Esto es una limitación en el sentido de que no podremos incluir esta variable (si existe o no netflix o hbo y el número de usuarios). Sin embargo sí que podremos acceder, no sin dificultad, a los datos del precio medio de una entrada del cine en España, en este link. En el gráfico superior podemos ver algunas estacionalidades, ya que los picos de mayor uso de las salas de cine son más o menos cíclicos. Lo comprobaremos a lo largo de la práctica observando las funciones de autocorrelación total y parcial, entre otros métodos que iremos viendo en el análisis.

Separaremos la muestra en datos de train y test, dejando como entrenamiento todos los datos previos al año 2009 y los 3 años restantes como período de prueba del modelo.

## Warning: package 'tseries' was built under R version 3.5.2

En los gráficos superiores podemos observar la descomposición de la serie en los diferentes elementos que la conforman. La tendencia en el tiempo es claramente descendente, mostrando la progresivamente menor popularidad de este medio de ocio. Podemos observar, además, que existe una cierta estacionalidad en los datos, pues los mismos “picos” de usuarios se repiten a lo largo del tiempo; el más grande de ellos se repite cada año. Es doble, pues hay al rededor de dos picos por año, con una bajada también grande antes de estos dos picos sucesivos. Respecto al ruido, parece oscilar al rededor de una media, aunque no es del todo constante en el tiempo en varianza, ya que hacia el año 2008 la varianza comienza a aumentar. Esto es positivo a la hora de aplicar un modelo GARCH, pues cuando las series tienen una homocedasticidad marcada este tipo de modelos no tienen tanto sentido, ya que la varianza del mes siguiente va a ser prácticamente igual a la del mes actual.

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(MASS)
## Warning: package 'MASS' was built under R version 3.5.2
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
box_cox <- boxcox(users ~ date,
                  data = train,
                  lambda = c(0, 0.5, 1))

lambda <- box_cox$x[which.max(box_cox$y)]
lambda
## [1] 0.6666667
powerTransform(train.ts)
## Estimated transformation parameter 
##  train.ts 
## 0.6566137

Vemos que \(\lambda\) está entre 0.5 y 1, por lo que tanto diferenciar como hacer la raíz cuadrada podrían ser buenas ideas. Lo que haremos más adelante será probar todo y quedarnos con el modelo que nos dé un menor error en la predicción, centrándonos más en el performance predictivo que en la robustness estadística del modelo, pues para ponerlo en producción lo que más dinero ahorrará o generará para la empresa será predecir de la mejor manera posible.

Funciones de Autocorrelación Parcial y Total: Ajuste de p y q.

Recordemos que la serie es mensual, por lo tanto el 1 del retardo hace referencia a un retardo de 12 meses (el eje X debemos multiplicarlo x 12 si queremos tenerlo en meses). Vemos en el gráfico de correlación parcial de la serie que existe una correlación significativa con el mes anterior, con 6 meses atrás, con 11 meses atrás y con 12 meses atrás. Esto haría referencia a la q del modelo, es decir la parte de la media móvil. Por otro lado vemos en la ACF que existe correlación con los datos de 1 año, 2 años y 3 años antes. También parece haber correlación con 1 mes antes, 8, 11 y 12 meses antes. Como vemos este comportamiento de que se correlacionan las observaciones con su homólogo de 1 año, 2 y 3 años antes, vamos a incluir una parte estacional.

Probamos primero a ajustar la parte estacional del modelo, de orden 1.

Primer Ajuste

## Series: train.ts 
## ARIMA(0,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##         sar1        mean
##       0.7340  10760.4226
## s.e.  0.0713    452.7668
## 
## sigma^2 estimated as 2397087:  log likelihood=-844.96
## AIC=1695.92   AICc=1696.18   BIC=1703.61
##                  sar1   intercept
## sar1       1.00000000 -0.04958451
## intercept -0.04958451  1.00000000
##      Retard    p-value
## [1,]      6 0.10506178
## [2,]     12 0.05530046
## [3,]     18 0.07325149
## [4,]     24 0.13259838
## [5,]     30 0.06617838
## [6,]     36 0.04317806
## [7,]     42 0.10556828
## [8,]     48 0.08014467

Vemos que se nos salen aún algunas correlaciones parciales de las bandas (de Barlett); nos llama especialmente en el PACF la banda del 1, ya que esta se sale mucho. En el caso de la ACF también se sale bastante esta misma banda. Como vemos en la PACF que existe correlación significativa con el valor de la serie en el mes anterior, decidimos probar un ARMA(0, 0, 1)x(0, 0, 1). Algo positivo es que el p-value de las correlaciones nos sale significativo para casi todos los lags, con lo cual estamos más cerca de cumplir con los supuestos del modelo ARMA.

Segundo Ajuste

ajuste2 <- Arima(train.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,0), period = 12),
                 method = "ML")

ajuste2
## Series: train.ts 
## ARIMA(0,0,1)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ma1    sar1        mean
##       0.3080  0.6986  10761.9363
## s.e.  0.1008  0.0770    518.2067
## 
## sigma^2 estimated as 2233367:  log likelihood=-840.47
## AIC=1688.95   AICc=1689.39   BIC=1699.21
##                    ma1        sar1    intercept
## ma1        1.000000000 -0.20087860  0.002181477
## sar1      -0.200878601  1.00000000 -0.047051282
## intercept  0.002181477 -0.04705128  1.000000000
##      Retard   p-value
## [1,]      6 0.7389242
## [2,]     12 0.2912638
## [3,]     18 0.5169280
## [4,]     24 0.7311593
## [5,]     30 0.4436429
## [6,]     36 0.3183765
## [7,]     42 0.4317232
## [8,]     48 0.4415686

Vemos que, especialmente en la PACF, persiste la correlación significativa con el valor del año anterior, por lo tanto decidimos incluir un nuevo componente estacional, y probaremos un ARMA(0, 0, 1)x(1, 0, 1).

Tercer Ajuste

ajuste3 <- Arima(train.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 method = "ML")

ajuste3
## Series: train.ts 
## ARIMA(0,0,1)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ma1    sar1     sma1        mean
##       0.3366  0.9012  -0.4260  10748.0593
## s.e.  0.1029  0.0492   0.1183    621.2377
## 
## sigma^2 estimated as 1984567:  log likelihood=-835.7
## AIC=1681.4   AICc=1682.07   BIC=1694.22
##                   ma1        sar1        sma1   intercept
## ma1        1.00000000  0.08343816 -0.25552459  0.00475227
## sar1       0.08343816  1.00000000 -0.68785287 -0.00188364
## sma1      -0.25552459 -0.68785287  1.00000000 -0.02603497
## intercept  0.00475227 -0.00188364 -0.02603497  1.00000000
##      Retard   p-value
## [1,]      6 0.2014748
## [2,]     12 0.1403626
## [3,]     18 0.2180477
## [4,]     24 0.4550048
## [5,]     30 0.3411850
## [6,]     36 0.2758098
## [7,]     42 0.4415517
## [8,]     48 0.4268327

Vemos que en la PACF se sigue saliendo una barra de los límites, en concreto la correlación de la observación en el período t con el valor de 2 años y 6 meses atrás… Podríamos probar algunos cambios en los órdenes del modelo para ver si esto desaparece, aunque no hay ninguno a primera vista que pudiera ser de ayuda. Una idea sería meterle una diferencia estacional. Una buena noticia es que seguimos sin tener que diferencir la serie obligatoriamente, pues los intervalos de confianza de nuestros coeficientes no incluyen al 1.

Cuarto Ajuste

Después de realizar muchas pruebas, nos encontramos con que, después del ARMA anterior que mostramos, el siguiente mejor sería un arima: ARMA(0, 1, 1)x(0, 1, 1). Sin embargo este modelo no mejoraría esto que comentábamos de que sigue existiendo una correlación parcial con las observaciones de dos años y medio atrás. Más adelante meteremos efectos de calendario y similares, y esperamos que estos nos puedan ayudar a explicar este suceso.

ajuste4 <- Arima(train.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(0,1,1), period = 12),
                 method = "ML")

ajuste4
## Series: train.ts 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.9185  -0.7404
## s.e.   0.0630   0.1364
## 
## sigma^2 estimated as 1563745:  log likelihood=-714.6
## AIC=1435.2   AICc=1435.51   BIC=1442.46
##              ma1        sma1
## ma1   1.00000000 -0.09393753
## sma1 -0.09393753  1.00000000
##      Retard   p-value
## [1,]      6 0.5054048
## [2,]     12 0.6491118
## [3,]     18 0.4477164
## [4,]     24 0.6099739
## [5,]     30 0.6912063
## [6,]     36 0.7900751
## [7,]     42 0.8462380
## [8,]     48 0.9300662

Probando con la raíz cuadrada:

Primer Ajuste

ajuste1_sqrt <- Arima(train2.ts,
                 order = c(0,0,0),
                 seasonal = list(order = c(1,0,0), period = 12),
                 method = "ML")

ajuste1_sqrt
## Series: train2.ts 
## ARIMA(0,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##         sar1      mean
##       0.7353  103.0976
## s.e.  0.0718    2.2057
## 
## sigma^2 estimated as 56.43:  log likelihood=-333.46
## AIC=672.92   AICc=673.18   BIC=680.61
##                sar1 intercept
## sar1       1.000000 -0.061306
## intercept -0.061306  1.000000
##      Retard    p-value
## [1,]      6 0.11317365
## [2,]     12 0.07335951
## [3,]     18 0.09586430
## [4,]     24 0.16246055
## [5,]     30 0.09605351
## [6,]     36 0.06446409
## [7,]     42 0.15921812
## [8,]     48 0.11712114

Segundo Ajuste

## Series: train2.ts 
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1    sar1      mean
##       0.3102  0.7042  103.1031
## s.e.  0.0982  0.0758    2.8135
## 
## sigma^2 estimated as 52.25:  log likelihood=-328.74
## AIC=665.48   AICc=665.92   BIC=675.74
##                    ar1        sar1    intercept
## ar1        1.000000000 -0.13106316 -0.001174168
## sar1      -0.131063164  1.00000000 -0.051907019
## intercept -0.001174168 -0.05190702  1.000000000
##      Retard   p-value
## [1,]      6 0.8491386
## [2,]     12 0.3958939
## [3,]     18 0.6000988
## [4,]     24 0.7914190
## [5,]     30 0.4284576
## [6,]     36 0.3247987
## [7,]     42 0.4542195
## [8,]     48 0.4616925

Tercer Ajuste

## Series: train2.ts 
## ARIMA(1,0,0)(1,1,1)[12] 
## 
## Coefficients:
##          ar1     sar1     sma1
##       0.3471  -0.3230  -0.2515
## s.e.  0.1157   0.1798   0.1898
## 
## sigma^2 estimated as 45.24:  log likelihood=-279.82
## AIC=567.65   AICc=568.15   BIC=577.37
##             ar1       sar1       sma1
## ar1   1.0000000  0.1754495 -0.4071704
## sar1  0.1754495  1.0000000 -0.7505946
## sma1 -0.4071704 -0.7505946  1.0000000
##      Retard    p-value
## [1,]      6 0.11419890
## [2,]     12 0.13121151
## [3,]     18 0.09610893
## [4,]     24 0.25061455
## [5,]     30 0.29791174
## [6,]     36 0.45922481
## [7,]     42 0.63002688
## [8,]     48 0.59753012

Parece haber una pequeña anomalía, pues encontramos una correlación parcial con las observaciones de hace año y medio significativa; sin embargo todo el resto de “barras” (correlaciones) se encuentran dentro de las bandas de Barlett.

Elevando la Serie a Lambda (Approximate Box-Cox Transformation)

Primer Ajuste

ajuste1_lambda <- Arima(train3.ts,
                 order = c(0,0,0),
                 seasonal = list(order = c(0,1,0), period = 12),
                 method = "ML")

ajuste1_lambda
## Series: train3.ts 
## ARIMA(0,0,0)(0,1,0)[12] 
## 
## sigma^2 estimated as 1.111:  log likelihood=-123.62
## AIC=249.23   AICc=249.28   BIC=251.66

Segundo Ajuste

ajuste2_lambda <- Arima(train3.ts,
                 order = c(0,0,0),
                 seasonal = list(order = c(0,1,1), period = 12),
                 method = "ML")

ajuste2_lambda
## Series: train3.ts 
## ARIMA(0,0,0)(0,1,1)[12] 
## 
## Coefficients:
##          sma1
##       -0.3652
## s.e.   0.0941
## 
## sigma^2 estimated as 0.9496:  log likelihood=-117.37
## AIC=238.74   AICc=238.89   BIC=243.6

Tercer Ajuste

Después de probar con varios, el modelo que mejor se ajusta es el SARIMA(0, 1, 1)x(1, 1, 1)[12]. Los residuos de este modelo presentan aún algunas anomalías, como veremos en los gráficos ACF y PACF, y es que no conseguimos que todas las correlaciones sean no significativas (dentro de la línea de Barlett) con \(y^\lambda, \lambda = 0.666667\). Esto podría ser debido a outliers, lo cual examinaremos más adelante.

ajuste3_lambda <- Arima(train3.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(1,1,1), period = 12),
                 method = "ML")

ajuste3_lambda
## Series: train3.ts 
## ARIMA(0,1,1)(1,1,1)[12] 
## 
## Coefficients:
##           ma1     sar1     sma1
##       -0.9043  -0.3177  -0.5540
## s.e.   0.0554   0.1595   0.1745
## 
## sigma^2 estimated as 0.6321:  log likelihood=-103.11
## AIC=214.21   AICc=214.72   BIC=223.89

Análisis de Intervenciones: Variables Explicativas.

Ahora vamos a añadirle a los mejores modelos de cada tipo (de cada transformación que hemos probado) las variables exógenas. Éstas van a ser relativas a si cae en semana santa, al número de días festivos del mes, al precio del cine medio en España ese año, y la variable dt, que se puede ver definida en el código.

En el chunk de código inferior vamos a definir una función que nos de las variables explicativas. Aparte de los festivos y demás efectos de calendario comunes a casi todas las series temporales, le añadiremos algunas concretas de este caso. Hemos cogido como “marco” la función definida por Daniel Vélez, y le hemos realizado diversos cambios para incorporar, por ejemplo, el efecto del día 2 de Mayo. Nos gustaría haber incluido que el día fuera miércoles como efecto, pues estos días hay descuento en los cines, pero es desde el año 2014 así que nos olvidamos. Sí vamos a incluir el precio medio del cine en cada año, pues esto esperamos que tenga una relación directa con el número de espectadores en las salas de España. Sacar estos datos nos llevará algo de trabajo, pues el único estudio disponible es este y no tiene una tabla de la que podamos sacar los datos, sino que estos tienen que ser sacados a mano. De igual manera, los introducimos como variables explicativas.

Primer Modelo

xregs = as.matrix(expl_train[ , c("semanaSanta", "dt", "bisiesto", "festivo", "price")])

ajuste5 <- Arima(train.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = xregs,
                 method = "ML")

ajuste5
## Series: train.ts 
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ma1    sar1     sma1  intercept  semanaSanta         dt  bisiesto
##       0.4353  0.9492  -0.6054  13895.825     360.7003  -113.3265  333.2582
## s.e.  0.0846  0.0362   0.1313   1313.897      74.9046    29.8079  690.8709
##        festivo      price
##       909.0092  -737.4628
## s.e.  350.3033   208.0835
## 
## sigma^2 estimated as 1270283:  log likelihood=-812.58
## AIC=1645.16   AICc=1647.75   BIC=1670.8
summary(ajuste5)
## Series: train.ts 
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ma1    sar1     sma1  intercept  semanaSanta         dt  bisiesto
##       0.4353  0.9492  -0.6054  13895.825     360.7003  -113.3265  333.2582
## s.e.  0.0846  0.0362   0.1313   1313.897      74.9046    29.8079  690.8709
##        festivo      price
##       909.0092  -737.4628
## s.e.  350.3033   208.0835
## 
## sigma^2 estimated as 1270283:  log likelihood=-812.58
## AIC=1645.16   AICc=1647.75   BIC=1670.8
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set -97.86562 1072.937 849.0645 -2.007611 8.058349 0.650576
##                    ACF1
## Training set 0.06290485
coeftest(ajuste5)
## 
## z test of coefficients:
## 
##                Estimate  Std. Error z value  Pr(>|z|)    
## ma1          4.3529e-01  8.4565e-02  5.1474 2.641e-07 ***
## sar1         9.4918e-01  3.6179e-02 26.2354 < 2.2e-16 ***
## sma1        -6.0538e-01  1.3127e-01 -4.6118 3.992e-06 ***
## intercept    1.3896e+04  1.3139e+03 10.5760 < 2.2e-16 ***
## semanaSanta  3.6070e+02  7.4905e+01  4.8155 1.469e-06 ***
## dt          -1.1333e+02  2.9808e+01 -3.8019 0.0001436 ***
## bisiesto     3.3326e+02  6.9087e+02  0.4824 0.6295403    
## festivo      9.0901e+02  3.5030e+02  2.5949 0.0094613 ** 
## price       -7.3746e+02  2.0808e+02 -3.5441 0.0003940 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box.test.2(residuals(ajuste5),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.7827383
## [2,]     12 0.8363230
## [3,]     18 0.7038717
## [4,]     24 0.8862800
## [5,]     30 0.9049797
## [6,]     36 0.8893456
## [7,]     42 0.8556661
## [8,]     48 0.9064435
cor.arma(ajuste5)
##                     ma1        sar1         sma1   intercept  semanaSanta
## ma1          1.00000000  0.03887024 -0.094470090 -0.04136217 -0.029139960
## sar1         0.03887024  1.00000000 -0.825996340  0.31795618  0.014912022
## sma1        -0.09447009 -0.82599634  1.000000000 -0.41203364  0.001940848
## intercept   -0.04136217  0.31795618 -0.412033638  1.00000000 -0.031848228
## semanaSanta -0.02913996  0.01491202  0.001940848 -0.03184823  1.000000000
## dt           0.07465806 -0.02273704  0.038908550 -0.01923734  0.142908868
## bisiesto     0.07091865  0.04288069 -0.094600353  0.06675116  0.055498652
## festivo     -0.01469069  0.04987521 -0.068116890 -0.19803690  0.058844929
## price        0.05357171 -0.34904465  0.450341177 -0.88010265 -0.008396677
##                      dt    bisiesto     festivo        price
## ma1          0.07465806  0.07091865 -0.01469069  0.053571711
## sar1        -0.02273704  0.04288069  0.04987521 -0.349044653
## sma1         0.03890855 -0.09460035 -0.06811689  0.450341177
## intercept   -0.01923734  0.06675116 -0.19803690 -0.880102650
## semanaSanta  0.14290887  0.05549865  0.05884493 -0.008396677
## dt           1.00000000 -0.01162851 -0.03982757  0.024539631
## bisiesto    -0.01162851  1.00000000  0.07225523 -0.107366327
## festivo     -0.03982757  0.07225523  1.00000000 -0.053487506
## price        0.02453963 -0.10736633 -0.05348751  1.000000000

Tenemos alguna correlación algo alta, como entre el precio y el intercept (-0.88); de todas formas la dejaremos de momento. Respecto a las demás variables, parece que la única que no es estadísticamente significativa es la binaria bisiesto. Una mala noticia es que, como vemos, el AIC baja poco con respecto al mismo modelo ARMA pero sin regresores externos. Vamos a volver a dibujar el ACF y PACF de los residuos del modelo, para comprobar si hemos conseguido resolver esa pequeña anomalía que observábamos (la correlación significativa con t - 2.5 años) (t - 30 meses, pues la serie es mensual). Antes de hacer esto debemos eliminar la variable no significativa (bisiesto).

Parece que sí ha quedado resuelto el pequeño conflicto que existía. Ahora ya tenemos todas las correlaciones dentro de las barras de Barlett y por lo tanto, teóricamente hemos capturado la información existente en la serie temporal.

Segundo Modelo: Con la Raíz Cuadrada

Como nos da problemas al sacar la diferencia estacional, cambiamos ligeramente el modelo presentado anteriormente con la Raíz Cuadrada. Usaremos un SARIMAX(1,0,0)x(1,0,1).

ajusteSqrt = Arima(train2.ts,
                   order=c(1,0,0),
                   seasonal=list(order=c(1,0,1), period=12),
                   xreg=xregs,
                   method="ML")

ajusteSqrt
## Series: train2.ts 
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sar1     sma1  intercept  semanaSanta       dt  bisiesto
##       0.4832  0.9583  -0.6101   117.2187       1.9221  -0.5387    0.8696
## s.e.  0.0915  0.0287   0.1206     7.9647       0.3536   0.1298    3.2162
##       festivo    price
##        4.4400  -3.3974
## s.e.   1.7961   1.2340
## 
## sigma^2 estimated as 27.68:  log likelihood=-298.3
## AIC=616.59   AICc=619.18   BIC=642.24
summary(ajusteSqrt)
## Series: train2.ts 
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sar1     sma1  intercept  semanaSanta       dt  bisiesto
##       0.4832  0.9583  -0.6101   117.2187       1.9221  -0.5387    0.8696
## s.e.  0.0915  0.0287   0.1206     7.9647       0.3536   0.1298    3.2162
##       festivo    price
##        4.4400  -3.3974
## s.e.   1.7961   1.2340
## 
## sigma^2 estimated as 27.68:  log likelihood=-298.3
## AIC=616.59   AICc=619.18   BIC=642.24
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.4333323 5.008137 3.950826 -0.6640638 3.855063 0.6270734
##                     ACF1
## Training set 0.008130705

Vemos en este caso que también conseguimos capturar todas las barras de correlaciones dentro de las bandas de Barlett, mejorando con los regresores externos el modelo expuesto anteriormente.

Tercer Modelo

Queremos probar también el otro modelo Arima sin transformaciones que habíamos planteado antes (0,1,1)x(0,1,1) con las variables exógenas. Esa era la intención, pero nos encontramos con un error persistente de R en el que nos dice “non-finite value supplied by optim”. Curiosamente este error sólo desaparece cuando quitamos la diferencia estacional: (1, 1, 1)x(1, 0, 1). Sin embargo esa diferencia estacional, una vez hemos decidido añadir SAr1, es necesaria… Por eso nos quedamos de momento con ese modelo tal cual, sin incluir las variables exógenas. Por otro lado, en el Jupyter Notebook adjunto en la carpeta que contiene este archivo vamos a intentar resolver este problema que nos hemos encontrado con Python, tratando de ver si éste puede optimizar la función que le especificamos.

En primer lugar, debe decirse que Python también tuvo problemas para optimizar el modelo SARIMAX que proponíamos, pero está algo mejor preparado el paquete statsmodels y tiene un método que le permite cambiar la fórmula utilizada para optimizar, sin utilizar gradientes ni hessianos, lo cual estaba causando problemas por algún motivo. Sin embargo, el ajuste es bastante malo (ver las estadísticas del Jupyter Notebook, en ellas se puede apreciar que casi ningún coeficiente sale significativo), una muestra de que existe un problema matemático al juntar el modelo con esa diferencia estacional con los regresores externos; al no poder utilizar los mejores métodos de optimización que son los que vienen por defecto, la estimación de los parámetros estará realizada de una forma más defectuosa, y por lo tanto el modelo no es tan bueno. Una muestra de esto es que, como vemos en los gráficos superiores, existe una autocorrelación parcial muy significativa con un período muy anterior, de muchos años… Esto es raro y es un comportamiento de los residuos que no hemos observado en el resto de modelos.

Cuarto Modelo.

Probaremos ahora con el modelo de \(y^\lambda\). De nuevo, quitamos la diferencia estacional para poder realizar los cálculos, pues con ella R no es capaz de optimizar la función y por lo tanto no puede estimar los parámetros.

ajusteLambda= Arima(train3.ts,
                   order=c(1,1,1),
                   seasonal=list(order=c(1,0,1), period=12),
                   xreg=xregs,
                   method="ML")

ajusteLambda
## Series: train3.ts 
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1  semanaSanta       dt  bisiesto
##       0.3718  -0.9443  0.9765  -0.6683       0.2617  -0.0714    0.1229
## s.e.  0.1049   0.0311  0.0184   0.1129       0.0476   0.0172    0.4205
##       festivo   price
##        0.5842  0.1647
## s.e.   0.2670  0.2538
## 
## sigma^2 estimated as 0.4488:  log likelihood=-101.77
## AIC=223.55   AICc=226.17   BIC=249.09
summary(ajusteLambda)
## Series: train3.ts 
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1  semanaSanta       dt  bisiesto
##       0.3718  -0.9443  0.9765  -0.6683       0.2617  -0.0714    0.1229
## s.e.  0.1049   0.0311  0.0184   0.1129       0.0476   0.0172    0.4205
##       festivo   price
##        0.5842  0.1647
## s.e.   0.2670  0.2538
## 
## sigma^2 estimated as 0.4488:  log likelihood=-101.77
## AIC=223.55   AICc=226.17   BIC=249.09
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.09179785 0.6340547 0.5019316 -0.5206396 2.413343 0.5958316
##                     ACF1
## Training set 0.003080306

En este caso los residuos sí presentan aún correlaciones significativas estadísticamente, por lo que no hemos conseguido mejorar demasiado en este caso introduciendo variables exógenas. Como vemos que hay algunas de estas variables que no son significativas (como el precio, por ejemplo), en este caso vamos a darle sólo algunas columnas de la matriz de regresores externos, las que pertenecen a las variables significativas.

## Series: train3.ts 
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1  semanaSanta       dt  festivo
##       0.3785  -0.9489  0.9744  -0.6598       0.2611  -0.0716   0.5841
## s.e.  0.1052   0.0312  0.0194   0.1129       0.0476   0.0173   0.2627
## 
## sigma^2 estimated as 0.4439:  log likelihood=-102.04
## AIC=220.07   AICc=221.75   BIC=240.5
## Series: train3.ts 
## Regression with ARIMA(1,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1  semanaSanta       dt  festivo
##       0.3785  -0.9489  0.9744  -0.6598       0.2611  -0.0716   0.5841
## s.e.  0.1052   0.0312  0.0194   0.1129       0.0476   0.0173   0.2627
## 
## sigma^2 estimated as 0.4439:  log likelihood=-102.04
## AIC=220.07   AICc=221.75   BIC=240.5
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.08939405 0.6378955 0.5088941 -0.5121413 2.445117 0.6040966
##                     ACF1
## Training set 0.002436342

Hay algo que venimos observando en la mayoría de los modelos, y es que cuando no introducimos una diferencia estacional el modelo lo pide a gritos, pues \(SAR1 \pm 2 \cdot sd(SAR1) > 1\).

Una vez hemos quitado las variables exógenas que estaban introduciendo ruido parece que el modelo se ajusta mucho mejor, capturando casi todas las autocorrelaciones tanto simples como parciales dentro de las barras de Barlett.

Quinto Modelo: Ajuste Automático

Probaremos también el ajuste automático del paquete forecast, para compararlo con los modelos que acabamos de construir.

## Series: train.ts 
## Regression with ARIMA(1,0,0)(1,0,0)[12] errors 
## 
## Coefficients:
##          ar1    sar1  intercept  semanaSanta         dt   bisiesto
##       0.4537  0.6961  12449.662     378.4617  -111.4025  -287.0549
## s.e.  0.0927  0.0775   1882.091      72.3948    34.9729   702.3373
##        festivo      price
##       840.6330  -467.3436
## s.e.  306.3595   317.3659
## 
## sigma^2 estimated as 1551111:  log likelihood=-820.35
## AIC=1658.7   AICc=1660.79   BIC=1681.77
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -58.6853 1192.414 918.0954 -1.832455 8.756608 0.7034693
##                    ACF1
## Training set 0.03301758
## 
## z test of coefficients:
## 
##                Estimate  Std. Error z value  Pr(>|z|)    
## ar1          4.5372e-01  9.2732e-02  4.8928 9.942e-07 ***
## sar1         6.9607e-01  7.7462e-02  8.9860 < 2.2e-16 ***
## intercept    1.2450e+04  1.8821e+03  6.6148 3.720e-11 ***
## semanaSanta  3.7846e+02  7.2395e+01  5.2277 1.716e-07 ***
## dt          -1.1140e+02  3.4973e+01 -3.1854  0.001446 ** 
## bisiesto    -2.8705e+02  7.0234e+02 -0.4087  0.682750    
## festivo      8.4063e+02  3.0636e+02  2.7439  0.006071 ** 
## price       -4.6734e+02  3.1737e+02 -1.4726  0.140867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##                      ar1        sar1   intercept  semanaSanta          dt
## ar1          1.000000000  0.02869996 -0.17252556  0.034373306  0.03105312
## sar1         0.028699963  1.00000000 -0.22470453  0.015800651  0.06786113
## intercept   -0.172525564 -0.22470453  1.00000000 -0.032907182 -0.02962391
## semanaSanta  0.034373306  0.01580065 -0.03290718  1.000000000  0.16617791
## dt           0.031053119  0.06786113 -0.02962391  0.166177908  1.00000000
## bisiesto    -0.032525076 -0.13584247  0.05270572  0.038990950  0.06340133
## festivo     -0.002152932 -0.08497891 -0.11188607  0.048950588 -0.03295701
## price        0.178976587  0.23922635 -0.93846931  0.008815508  0.02856165
##                bisiesto      festivo        price
## ar1         -0.03252508 -0.002152932  0.178976587
## sar1        -0.13584247 -0.084978910  0.239226352
## intercept    0.05270572 -0.111886071 -0.938469314
## semanaSanta  0.03899095  0.048950588  0.008815508
## dt           0.06340133 -0.032957006  0.028561651
## bisiesto     1.00000000  0.105985726 -0.083067380
## festivo      0.10598573  1.000000000 -0.043296371
## price       -0.08306738 -0.043296371  1.000000000
##      Retard   p-value
## [1,]      6 0.9712270
## [2,]     12 0.5183626
## [3,]     18 0.6650078
## [4,]     24 0.6116338
## [5,]     30 0.5438551
## [6,]     36 0.5760761
## [7,]     42 0.6409670
## [8,]     48 0.6901314

En este caso el precio parece no ser significativo, tendríamos que coger \(\alpha = 0.15\) para que éste fuera significativo. Esto podría estar introduciendo algo de ruido… El Stepwise lo coge porque seguramente el AIC baje al introducirlo, pero no nos convence.

ajusteAutomatico2<-auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1, 
                             seasonal=TRUE,ic="aic",xreg=xregs[ , 1:4],stepwise=TRUE)

summary(ajusteAutomatico2)
## Series: train.ts 
## Regression with ARIMA(1,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1  semanaSanta         dt   bisiesto   festivo
##       0.4309  -0.9603  0.7183     378.2619  -109.6057  -306.8391  833.8064
## s.e.  0.0998   0.0254  0.0748      71.3012    34.5435   688.3324  319.6909
## 
## sigma^2 estimated as 1534906:  log likelihood=-812.37
## AIC=1640.74   AICc=1642.41   BIC=1661.17
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -155.7495 1186.169 917.4091 -2.598247 8.826862 0.7029434
##                    ACF1
## Training set 0.01743244
coeftest(ajusteAutomatico2)
## 
## z test of coefficients:
## 
##                Estimate  Std. Error  z value  Pr(>|z|)    
## ar1            0.430888    0.099837   4.3159 1.589e-05 ***
## ma1           -0.960282    0.025393 -37.8173 < 2.2e-16 ***
## sar1           0.718279    0.074777   9.6056 < 2.2e-16 ***
## semanaSanta  378.261909   71.301206   5.3051 1.126e-07 ***
## dt          -109.605727   34.543479  -3.1730  0.001509 ** 
## bisiesto    -306.839073  688.332416  -0.4458  0.655762    
## festivo      833.806369  319.690888   2.6082  0.009103 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajusteAutomatico2)
##                      ar1         ma1         sar1 semanaSanta          dt
## ar1          1.000000000 -0.32387322 -0.007077443  0.03629576  0.02554625
## ma1         -0.323873218  1.00000000 -0.202175178 -0.01596965 -0.01307942
## sar1        -0.007077443 -0.20217518  1.000000000  0.02190270  0.06614748
## semanaSanta  0.036295761 -0.01596965  0.021902700  1.00000000  0.16615284
## dt           0.025546254 -0.01307942  0.066147478  0.16615284  1.00000000
## bisiesto    -0.031243024  0.04894940 -0.115393743  0.03544888  0.06710558
## festivo     -0.003781187  0.02558829 -0.069071152  0.04473960 -0.03011573
##                bisiesto      festivo
## ar1         -0.03124302 -0.003781187
## ma1          0.04894940  0.025588287
## sar1        -0.11539374 -0.069071152
## semanaSanta  0.03544888  0.044739603
## dt           0.06710558 -0.030115726
## bisiesto     1.00000000  0.097608925
## festivo      0.09760893  1.000000000
Box.test.2(residuals(ajusteAutomatico2),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.9244382
## [2,]     12 0.2802904
## [3,]     18 0.4741342
## [4,]     24 0.4239821
## [5,]     30 0.3697063
## [6,]     36 0.4117417
## [7,]     42 0.4889451
## [8,]     48 0.5470928
acf(ajusteAutomatico2$residuals, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(ajusteAutomatico2$residuals, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

Vemos que no cambia demasiado con respecto a cuando metemos price como regresor externo. Nos quedamos con el primero, pues tiene un MAPE algo menor.

Identificación de Outliers

En esta sección nos vamos a encargar de detectar los outliers y las anomalías de los modelos que hemos empleado en R, pues el que hemos empleado en Python tendría que ser tratado con paquetes de Python y con otro tipo de procedimientos. Nos centramos por tanto en los siguientes modelos:

  • Primer Modelo: SARIMAX(0, 0, 1)x(1, 0, 1)[12]

  • Segundo Modelo: \(\sqrt{y}\) SARIMAX(1, 0, 0)x(1, 0, 1)[12]

  • Tercer Modelo: \(y^\lambda\) SARIMAX (0, 1, 1)x(1, 0, 1)[12]

  • Cuarto Modelo: Ajuste Automático SARIMAX(1,1,1)x(1,0,0)[12]

En este apartado utilizaremos siglas para referirnos a los diferentes tipos de anomalías o outliers que existen en una serie temporal, para abreviar el texto. Por ello, aquí vamos a especificar de qué tipo de outliers hablamos cuando utilizamos cada una de estas siglas.

  • AO: Additive Outliers

  • IO: Innovative Outliers: estos son un tipo de outliers que dependen del modelo que se esté utilizando, en contraste con el resto. Más información sobre este tema, así sobre el resto de tipos de outliers aquí.

  • LS: Level Shift.

  • TC: Transitory Change.

El propio nombre de cada tipo de outlier es suficientemente informativo, ya que nos dice cuál es el tipo de cambio que sufre la serie en esa observación que consideramos outlier.

Identificación de los Outliers del Primer Modelo

## Warning: package 'tsoutliers' was built under R version 3.5.2
##   ind type    coefhat     tstat abststat       date
## 1  43   TC  3221.6158  3.998365 3.998365 2004-07-01
## 2  67   LS  -992.9195 -3.040895 3.040895 2006-07-01
## 3  70   LS -1007.7600 -3.036758 3.036758 2006-10-01
## 4  81   AO -2837.9402 -3.356161 3.356161 2007-09-01

Vaya, parece que si encontramos algunos patrones de outliers. En concreto, vemos que tenemos un outlier repetido, aunque de diferente tipo (el primero es de tipo TC y el segundo de tipo LS) el 1 de Julio de 2004 y 2006. Por otro lado, el 1 de Octubre de 2006 tenemos otro del tipo LS y el 1 de Septiembre del 2007 tenemos uno del tipo AO.

library(ggplot2)

v = rep(0, nrow(train))

v[c(40:46, 66:68, 69:71, 78:84)] <- 1

inds = diff(c(0, v))
start = train$date[inds == 1]
end = train$date[inds==-1]

if(length(start) > length(end)) end <- c(end, tail(train$date, 1))

rects = data.frame(start=start, end = end, group=seq_along(start))

ggplot(data = train, aes(date, users)) +
  theme_minimal() +
  geom_line(lty=2, color="steelblue", lwd=1.1) +
  
  geom_point() +
  
  geom_rect(data = rects, inherit.aes=FALSE, aes(xmin=start, xmax=end, ymin=min(train$users),
                                                 ymax=max(train$users), group=group), color="transparent", fill="orange", alpha=0.3) +
  
    geom_vline(aes(xintercept=as.Date("2007-09-01"))) +

  geom_vline(aes(xintercept=as.Date("2004-07-01"))) +
  
    geom_vline(aes(xintercept=as.Date("2006-07-01"))) +

    geom_vline(aes(xintercept=as.Date("2006-10-01"))) +

  ggtitle("Outliers of the Train Set with the First Model")

En el gráfico superior podemos ver los diferentes outliers que presenta el primer modelo.

##   ind type    coefhat     tstat abststat       date
## 1  43   TC  3221.6158  3.998365 3.998365 2004-07-01
## 2  81   AO -2837.9402 -3.356161 3.356161 2007-09-01
## 3  67   LS  -992.9195 -3.040895 3.040895 2006-07-01
## 4  70   LS -1007.7600 -3.036758 3.036758 2006-10-01

Ahora ajustaremos de nuevo el modelo teniendo en cuenta los outliers. Vemos que las variables LS67 y LS70 no son estadísticamente significativas, por lo que las eliminamos. Eliminamos también bisiesto que no es significativa.

ajuste5def = Arima(train.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = RegressorsMasOutliers[, c(1, 2, 4,5,6,7)],
                 method = "ML")

ajuste5def
## Series: train.ts 
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ma1    sar1     sma1  intercept  semanaSanta         dt   festivo
##       0.4718  0.9553  -0.5802  13902.622     356.4204  -127.0938  862.4788
## s.e.  0.0929  0.0296   0.1263   1197.981      61.6411    26.0868  335.2374
##           price       TC43        AO81
##       -737.9075  3133.9434  -2871.3002
## s.e.   186.0039   757.3415    786.1132
## 
## sigma^2 estimated as 951818:  log likelihood=-799.23
## AIC=1620.45   AICc=1623.6   BIC=1648.66
summary(ajuste5def)
## Series: train.ts 
## Regression with ARIMA(0,0,1)(1,0,1)[12] errors 
## 
## Coefficients:
##          ma1    sar1     sma1  intercept  semanaSanta         dt   festivo
##       0.4718  0.9553  -0.5802  13902.622     356.4204  -127.0938  862.4788
## s.e.  0.0929  0.0296   0.1263   1197.981      61.6411    26.0868  335.2374
##           price       TC43        AO81
##       -737.9075  3133.9434  -2871.3002
## s.e.   186.0039   757.3415    786.1132
## 
## sigma^2 estimated as 951818:  log likelihood=-799.23
## AIC=1620.45   AICc=1623.6   BIC=1648.66
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -94.48002 923.4014 726.6795 -1.659199 6.872925 0.5568013
##                    ACF1
## Training set 0.02273563

Identificación de Outliers del Segundo Modelo

##   ind type   coefhat     tstat abststat       date
## 1  81   AO -14.58462 -3.415948 3.415948 2007-09-01

En este caso solamente tenemos un outlier del tipo AO el 1 de Septiembre de 2007.

En este caso eliminamos la variable bisiesto pues no es significativa, y ajustamos el siguiente modelo:

## Series: train2.ts 
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sar1     sma1  intercept  semanaSanta       dt  festivo
##       0.5230  0.9618  -0.5908   116.1457       1.8731  -0.6526   4.2460
## s.e.  0.0901  0.0260   0.1241     8.2292       0.3164   0.1201   1.7845
##         price      AO81
##       -3.1451  -15.4789
## s.e.   1.2509    3.9725
## 
## sigma^2 estimated as 23.6:  log likelihood=-291.45
## AIC=602.89   AICc=605.48   BIC=628.54
## Series: train2.ts 
## Regression with ARIMA(1,0,0)(1,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sar1     sma1  intercept  semanaSanta       dt  festivo
##       0.5230  0.9618  -0.5908   116.1457       1.8731  -0.6526   4.2460
## s.e.  0.0901  0.0260   0.1241     8.2292       0.3164   0.1201   1.7845
##         price      AO81
##       -3.1451  -15.4789
## s.e.   1.2509    3.9725
## 
## sigma^2 estimated as 23.6:  log likelihood=-291.45
## AIC=602.89   AICc=605.48   BIC=628.54
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.4389988 4.625089 3.669041 -0.6284198 3.555108 0.5823486
##                      ACF1
## Training set -0.004087827
## 
## z test of coefficients:
## 
##               Estimate Std. Error z value  Pr(>|z|)    
## ar1           0.523007   0.090058  5.8075 6.343e-09 ***
## sar1          0.961818   0.026033 36.9464 < 2.2e-16 ***
## sma1         -0.590811   0.124114 -4.7602 1.934e-06 ***
## intercept   116.145686   8.229197 14.1139 < 2.2e-16 ***
## semanaSanta   1.873100   0.316407  5.9199 3.221e-09 ***
## dt           -0.652622   0.120132 -5.4325 5.556e-08 ***
## festivo       4.245998   1.784493  2.3794   0.01734 *  
## price        -3.145081   1.250873 -2.5143   0.01193 *  
## AO81        -15.478855   3.972484 -3.8965 9.759e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Identificación de Outliers del Tercer Modelo

##   ind type   coefhat     tstat abststat       date
## 1  81   AO -1.924749 -3.610469 3.610469 2007-09-01

Sale exactamente el mismo outlier que en el caso anterior (también había aparecido en el primero). Esto contradice lo que se mencionaba en el link que se encuentra más arriba, en el que se decía que los outliers del tipo IO no dependen del modelo a diferencia de los demás. Por lo que estamos observando, el outlier que no depende del modelo es el AO de 2007-09-01.

## $title
## [1] "Outliers of the Train Set with the Third Model"
## 
## attr(,"class")
## [1] "labels"

En este caso eliminamos bisiesto y price, pues no son significativas.

ajusteLambdaDef = Arima(train3.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = RegressorsMasOutliers3[, c(1, 2, 4, 6)],
                 method = "ML")

ajusteLambdaDef
## Series: train3.ts 
## Regression with ARIMA(0,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##           ma1    sar1     sma1  semanaSanta       dt  festivo     AO81
##       -0.5740  0.9820  -0.6699       0.2548  -0.0900   0.5410  -2.2398
## s.e.   0.1552  0.0147   0.1198       0.0478   0.0178   0.2848   0.5493
## 
## sigma^2 estimated as 0.4275:  log likelihood=-101.68
## AIC=219.36   AICc=221.03   BIC=239.79
summary(ajusteLambdaDef)
## Series: train3.ts 
## Regression with ARIMA(0,1,1)(1,0,1)[12] errors 
## 
## Coefficients:
##           ma1    sar1     sma1  semanaSanta       dt  festivo     AO81
##       -0.5740  0.9820  -0.6699       0.2548  -0.0900   0.5410  -2.2398
## s.e.   0.1552  0.0147   0.1198       0.0478   0.0178   0.2848   0.5493
## 
## sigma^2 estimated as 0.4275:  log likelihood=-101.68
## AIC=219.36   AICc=221.03   BIC=239.79
## 
## Training set error measures:
##                       ME      RMSE      MAE        MPE    MAPE      MASE
## Training set -0.02216439 0.6260056 0.508905 -0.1635231 2.43026 0.6041096
##                   ACF1
## Training set 0.1444174
coeftest(ajusteLambdaDef)
## 
## z test of coefficients:
## 
##              Estimate Std. Error z value  Pr(>|z|)    
## ma1         -0.574004   0.155191 -3.6987 0.0002167 ***
## sar1         0.981959   0.014662 66.9715 < 2.2e-16 ***
## sma1        -0.669914   0.119799 -5.5920 2.245e-08 ***
## semanaSanta  0.254823   0.047787  5.3324 9.691e-08 ***
## dt          -0.089971   0.017800 -5.0545 4.315e-07 ***
## festivo      0.541018   0.284765  1.8999 0.0574490 .  
## AO81        -2.239846   0.549277 -4.0778 4.546e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Identificación de Outliers del Cuarto Modelo

##   ind type   coefhat     tstat abststat       date
## 1  26   AO  2395.624  3.290378 3.290378 2003-02-01
## 2  40   TC  2295.855  3.029250 3.029250 2004-04-01
## 3  43   TC  2959.651  3.905092 3.905092 2004-07-01
## 4  81   AO -2919.068 -4.009327 4.009327 2007-09-01

En este caso observamos unos outliers algo distintos, en concreto el del 1 de Frebrero del 2003.

## Series: train.ts 
## Regression with ARIMA(0,0,1)(1,0,0)[12] errors 
## 
## Coefficients:
##          ma1    sar1  intercept  semanaSanta         dt  bisiesto
##       0.5725  0.7442  12839.046     323.7978  -144.1346  985.1539
## s.e.  0.1154  0.0742   1398.947      50.7277    29.1279  564.1778
##        festivo      price       AO26       TC43        AO81
##       862.9701  -555.2582  2915.7790  3174.1653  -3092.7717
## s.e.  257.3901   234.4697   710.8307   718.1577    664.0621
## 
## sigma^2 estimated as 1047564:  log likelihood=-800.79
## AIC=1625.58   AICc=1629.34   BIC=1656.35
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set -67.8566 963.0838 778.422 -1.539647 7.378763 0.5964477
##                     ACF1
## Training set 0.000512715
## 
## z test of coefficients:
## 
##                Estimate  Std. Error z value  Pr(>|z|)    
## ma1          5.7247e-01  1.1541e-01  4.9602 7.041e-07 ***
## sar1         7.4415e-01  7.4153e-02 10.0354 < 2.2e-16 ***
## intercept    1.2839e+04  1.3989e+03  9.1777 < 2.2e-16 ***
## semanaSanta  3.2380e+02  5.0728e+01  6.3831 1.736e-10 ***
## dt          -1.4413e+02  2.9128e+01 -4.9483 7.485e-07 ***
## bisiesto     9.8515e+02  5.6418e+02  1.7462 0.0807803 .  
## festivo      8.6297e+02  2.5739e+02  3.3528 0.0008001 ***
## price       -5.5526e+02  2.3447e+02 -2.3681 0.0178775 *  
## AO26         2.9158e+03  7.1083e+02  4.1019 4.097e-05 ***
## TC43         3.1742e+03  7.1816e+02  4.4199 9.876e-06 ***
## AO81        -3.0928e+03  6.6406e+02 -4.6574 3.203e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##                      ma1        sar1    intercept  semanaSanta
## ma1          1.000000000 -0.20437116 -0.135490559 -0.203795824
## sar1        -0.204371158  1.00000000 -0.206725827  0.068018944
## intercept   -0.135490559 -0.20672583  1.000000000 -0.013153707
## semanaSanta -0.203795824  0.06801894 -0.013153707  1.000000000
## dt          -0.005670653  0.06776805 -0.045267277  0.203753673
## bisiesto     0.261945821 -0.15536832  0.048395787 -0.091535287
## festivo     -0.018197633 -0.09188222 -0.124342996  0.032923276
## price        0.146219717  0.22662111 -0.928634546 -0.007164621
## AO26         0.095239154 -0.09570366  0.123411439 -0.174263687
## TC43         0.220047765 -0.12254738 -0.001739071  0.006132956
## AO81        -0.068885427  0.04353456 -0.054292399  0.077101984
##                       dt    bisiesto      festivo        price
## ma1         -0.005670653  0.26194582 -0.018197633  0.146219717
## sar1         0.067768054 -0.15536832 -0.091882220  0.226621113
## intercept   -0.045267277  0.04839579 -0.124342996 -0.928634546
## semanaSanta  0.203753673 -0.09153529  0.032923276 -0.007164621
## dt           1.000000000 -0.08464818 -0.067611369  0.052000771
## bisiesto    -0.084648181  1.00000000  0.111077678 -0.083761033
## festivo     -0.067611369  0.11107768  1.000000000 -0.050783867
## price        0.052000771 -0.08376103 -0.050783867  1.000000000
## AO26        -0.146432322  0.40544776  0.070151245 -0.148718697
## TC43         0.060416959  0.09769199  0.006446041 -0.010730443
## AO81         0.277853894 -0.02177463 -0.002048460  0.051944791
##                     AO26         TC43         AO81
## ma1          0.095239154  0.220047765 -0.068885427
## sar1        -0.095703658 -0.122547377  0.043534558
## intercept    0.123411439 -0.001739071 -0.054292399
## semanaSanta -0.174263687  0.006132956  0.077101984
## dt          -0.146432322  0.060416959  0.277853894
## bisiesto     0.405447761  0.097691991 -0.021774634
## festivo      0.070151245  0.006446041 -0.002048460
## price       -0.148718697 -0.010730443  0.051944791
## AO26         1.000000000  0.009820674 -0.045422480
## TC43         0.009820674  1.000000000  0.001388494
## AO81        -0.045422480  0.001388494  1.000000000
##      Retard    p-value
## [1,]      6 0.10212090
## [2,]     12 0.03233052
## [3,]     18 0.01566603
## [4,]     24 0.00889139
## [5,]     30 0.00779644
## [6,]     36 0.01289081
## [7,]     42 0.01397820
## [8,]     48 0.02218604

Vemos que la serie sigue presentando anomalías. Además, el ajuste automático ha cambiado, y ahora estima un SARIMAX(0, 0, 1)x(1, 0, 0).

Análisis de la Varianza

Primer Modelo

La pregunta que nos hacemos antes de ponernos a predecir es la siguiente, ¿admitiría alguno de los modelos que hemos visto arriba el ajuste de un modelo GARCH para modelizar la varianza de la serie?.

##      Retard   p-value
## [1,]      6 0.6938020
## [2,]     12 0.8881881
## [3,]     18 0.7590623
## [4,]     24 0.8721452
## [5,]     30 0.9708096
## [6,]     36 0.9898367
## [7,]     42 0.9969951
## [8,]     48 0.9829229

En este caso los residuos presentan ruido blanco, por lo tanto no vamos a poder utilizar un modelo GARCH.

Segundo Modelo

##      Retard    p-value
## [1,]      6 0.05882292
## [2,]     12 0.12071514
## [3,]     18 0.02685252
## [4,]     24 0.07757110
## [5,]     30 0.07584400
## [6,]     36 0.17743292
## [7,]     42 0.21977585
## [8,]     48 0.28055270

En este caso el ruido blanco no está tan claro, por lo que sí podría tener sentido utilizar un modelo GARCH para modelizar la varianza de la serie. Sin embargo, la modelización de la varianza de la serie temporal va a ser llevada acabo únicamente con la serie ajustada de manera automática, pues lo vamos a utilizar para contrastar un “Machine Learning Approach” versus los ajustes manuales que hace el humano con su conocimiento de estadística y su detección de patrones (al final lo que hacemos al analizar los gráficos ACF y PACF y ajustando de nuevo los parámetros es detectar patrones en la serie temporal y jugar con ellos).

Tercer Modelo

##      Retard    p-value
## [1,]      6 0.02051878
## [2,]     12 0.09414880
## [3,]     18 0.26736440
## [4,]     24 0.30652493
## [5,]     30 0.21065611
## [6,]     36 0.31032755
## [7,]     42 0.52150333
## [8,]     48 0.58652776

Algo menos que en el anterior caso, quizás sí se podría utilizar un GARCH.

Cuarto Modelo

##      Retard   p-value
## [1,]      6 0.3452732
## [2,]     12 0.5060421
## [3,]     18 0.3245441
## [4,]     24 0.4997114
## [5,]     30 0.7085126
## [6,]     36 0.8029283
## [7,]     42 0.2567279
## [8,]     48 0.3872969

En este caso queda claro que el proceso es de ruido blanco y por lo tanto no parece tener mucho sentido ajustar un modelo GARCH. Lo que vamos a hacer en la siguiente sección será ajustar primero los 3 modelos manuales únicamente para el modelo de la media, y por otro lado, con una función diferente, vamos a realizar un ajuste automático en el que simultáneamente predigamos el valor medio, así como los valores superiores e inferiores, obtenidos estos mediante la predicción de la desviación estándar con un modelo GARCH(1,1).

PREDICCIONES, EVALUACIÓN DE MODELOS Y RESULTADOS

Time Series Approach

get_predictions = function(d = datos, nendtrain=96, var='users', datecol='date', lam=0.6566137){
  
  ########################################################################################################################
  #     This function will calculate the predictions, in periods of 12 months, for the cinema users in Spain            #
  #     with the different models presented.                                                                            #
  #     Parameters:                                                                                                     #
  #       - d: data                                                                                                     #
  #       - nendtrain: index of the last observation for the pure train part, before starting to predict (moving window)#
  #       - var: the variable of interest inside the dataset for univariate time series                                 #
  #       - datecol: the name or index of the date variable in the dataset                                              #     
  #       - lam: the lambda value for Box-Cox transformation as returned by car::powerTransform()                       #
  #     Returns:                                                                                                        #
  #       - preds_df: data frame with predictions and real value.                                                       #
  #######################################################################################################################
  
  preds1 = c()
  
  preds2 = c()
  
  preds3 = c()
  
  for (i in seq(nendtrain, nrow(d) - 12, 12)){
    
    #dates = d[1:i, datecol]
    
    ###########################
    #     TRAIN & TEST        #             
    #       SPLITS            #
    ###########################
    
    tr1 = data.frame(d[1:i,])
    
    tr2 = data.frame(d[1:i, ])
    
    tr2[, var] = sqrt(tr2[, var])
    
    tr3 = data.frame(d[1:i, ])
    
    tr3[, var] = tr3[,var]**lam
    
    te1 = data.frame(d[(i+1):(i+12), ])
    
    te2 = data.frame(d[(i+1):(i+12), ])
    
    te2[,var] = sqrt(te2[,var])
    
    te3 = data.frame(d[(i+1):(i+12), ])
    
    te3[,var] = te3[,var]**lam
    
    #####################
    #     COMPUTING     #
    # EXTERNAL VARIABLES#
    #####################
    
    xregs = as.matrix(calculoExplicativasCalendario(tr1[,datecol],0)[ , c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
    
    newxregs = as.matrix(calculoExplicativasCalendario(te1$date, 0)[,c("semanaSanta","dt","bisiesto","festivo","price")])

    tr1.ts = ts(tr1[,var], start=c(year(tr1[1, datecol]), month(tr1[1, datecol])), 
                   end=c(year(tr1[nrow(tr1), datecol]), month(tr1[nrow(tr1), datecol])),
                   frequency=12)
    
    tr2.ts = ts(tr2[,var], start=c(year(tr2[1, datecol]), month(tr2[1, datecol])), 
                   end=c(year(tr2[nrow(tr2), datecol]), month(tr2[nrow(tr2), datecol])),
                   frequency=12)
    
    tr3.ts = ts(tr3[,var], start=c(year(tr3[1, datecol]), month(tr3[1, datecol])), 
                   end=c(year(tr3[nrow(tr3), datecol]), month(tr3[nrow(tr3), datecol])),
                   frequency=12)
    
    te1.ts = ts(te1[, var], start=c(year(te1[1, datecol]), month(te1[1,datecol])),
                  end = c(year(te1[nrow(te1), datecol]), month(te1[nrow(te1),datecol])), 
                  frequency=12)
    
    te2.ts = ts(te2[, var], start=c(year(te2[1, datecol]), month(te2[1,datecol])),
                  end = c(year(te2[nrow(te2), datecol]), month(te2[nrow(te2),datecol])), 
                  frequency=12)
    
    te3.ts = ts(te3[, var], start=c(year(te3[1, datecol]), month(te3[1,datecol])),
                  end = c(year(te3[nrow(te3), datecol]), month(te3[nrow(te3),datecol])), 
                  frequency=12)
    
    ###############################################################################################################
    #   FITTING, OUTLIERS CALCULATION, EXTERNAL VARIABLES RECALCULATED WITH OUTLIERS; PREDICTION FOR EACH MODEL   #
    ###############################################################################################################
    
    
    #################
    #   1ST MODEL   #
    #################
    
    ajuste1 <- Arima(tr1.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = xregs,
                 method = "ML")
    
    listaOutliers1 <- locate.outliers(ajuste1$residuals,
                                 pars = coefs2poly(ajuste1),
                                 types = c("AO", "LS", "TC"),cval=3)

   
    outliers1 <- outliers(c(as.character(listaOutliers1$type)), c(listaOutliers1$ind))
    outliersVariables1 <- outliers.effects(outliers1, length(ajuste1$residuals))
    RegressorsMasOutliers1 <- cbind(xregs,outliersVariables1)
    
    ajuste1def=Arima(tr1.ts,
                 order = c(0,0,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = RegressorsMasOutliers1[, -3],
                 method = "ML")
    
    
    outliersNew1 <- outliersVariables1[(i-11):i,]

    newxregs1 <- as.matrix(cbind(newxregs,outliersNew1))

    prediccion1 <- predict(ajuste1def, n.ahead=12,
                          newxreg=newxregs1[, -3])
    
    preds1 = c(preds1, prediccion1$pred)
    
    #################
    #   2ND MODEL   #
    #################
    
    ajuste2 = Arima(tr2.ts,
                   order=c(1,0,0),
                   seasonal=list(order=c(1,0,1), period=12),
                   xreg=xregs,
                   method="ML")
    
    listaOutliers2 <- locate.outliers(ajuste2$residuals,
                                 pars = coefs2poly(ajuste2),
                                 types = c("AO", "LS", "TC"),cval=3)
    
    outliers2 <- outliers(c(as.character(listaOutliers2$type)), c(listaOutliers2$ind))
    
    tryCatch(
      {outliersVariables2 <- outliers.effects(outliers2, length(ajuste2$residuals))
    RegressorsMasOutliers2 <- cbind(xregs,outliersVariables2)
    },
    error=function(e){RegressorsMasOutliers2 = as.matrix(xregs)}
    )
    
    if(nrow(RegressorsMasOutliers2 != length(tr2.ts))){
      print("check this problem again")
      a=0
      RegressorsMasOutliers2=as.matrix(xregs)
    }else{a = 1}
    
    ajuste2def = Arima(tr2.ts,
                 order = c(1,0,0),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = RegressorsMasOutliers2[, -3],
                 method = "ML")
    
    if(a == 0){
      
      newxregs2 <- as.matrix(newxregs[,-3])
      
      prediccion2 <- predict(ajuste2def, n.ahead=12,
                             newxreg=newxregs2)
      
      preds2 = c(preds2, prediccion2$pred)
      
    }else{
      
    outliersNew2 <- outliersVariables2[(i-11):i,]

    newxregs2 <- as.matrix(cbind(newxregs,outliersNew2))

    prediccion2 <- predict(ajuste2def, n.ahead=12,
                          newxreg=newxregs2[, -3])
    
    preds2 = c(preds2, prediccion2$pred)
      
    }
    

    
    ###############
    #   3RD MODEL #
    ###############
    
    ajuste3 = Arima(tr3.ts,
                   order=c(1,1,1),
                   seasonal=list(order=c(1,0,1), period=12),
                   xreg=xregs[, -c(3,5)],
                   method="ML")
    
   listaOutliers3 <- locate.outliers(ajuste3$residuals,
                                 pars = coefs2poly(ajuste3),
                                 types = c("AO", "LS", "TC"),cval=3)
    
    outliers3 <- outliers(c(as.character(listaOutliers3$type)), c(listaOutliers3$ind))
    outliersVariables3 <- outliers.effects(outliers3, length(ajuste3$residuals))
    RegressorsMasOutliers3 <- cbind(xregs,outliersVariables3)
    
    ajuste3def=Arima(tr3.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(1,0,1), period = 12),
                 xreg = RegressorsMasOutliers3[, -c(3,5)],
                 method = "ML")
    
    outliersNew3 <- outliersVariables3[(i-11):i,]

    newxregs3 <- as.matrix(cbind(newxregs,outliersNew3))

    prediccion3 <- predict(ajuste3def, n.ahead=12,
                          newxreg=newxregs3[, -c(3,5)])
    
    preds3 = c(preds3, prediccion3$pred) 
    
  }
  
  if(length(preds1)!=length(preds2) | length(preds1) != length(preds3)){
    return(list(preds1, preds2, preds3))
  }
  
  real = data.frame(d)[(nendtrain+1):nrow(d), var]
  
  preds_df = data.frame(modelo1=preds1,modelo2=preds2,modelo3=preds3, real=real)
  
  return(preds_df)
  
}
## [1] "check this problem again"
## [1] "check this problem again"
## [1] "check this problem again"
## [1] "check this problem again"
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Como podemos ver, los modelos 1 y 2 se “pisan” prácticamente, clavando el uno las predicciones del otro. Estos dos tienden a estimar algo por encima del valor real, como vemos en el gráfico de arriba, mientras que el tercer modelo suele estimar algo por debajo. Visualmente parece que el modelo que mejor predice es el 3; vamos a constatarlo en números.

Estadísticas de Error del Modelo 1

##                 ME     RMSE      MAE       MPE     MAPE
## Test set -809.9391 1527.874 1208.677 -11.54855 15.44051

Estadísticas de Error del Modelo 2

##                 ME     RMSE      MAE       MPE     MAPE
## Test set -790.9721 1463.625 1153.604 -11.20356 14.70737

Estadísticas de Error del Modelo 3

##                ME    RMSE      MAE      MPE     MAPE
## Test set 297.8611 1222.87 996.9501 2.331691 11.63035

Como decíamos, el modelo 3 es el que mejor predice de estos 3, ya que tiene un enor MAPE y menor RMSE (el MAPE nos permite medir el error en términos porcentuales, en este caso un 11.63% (está multiplicado por 100)); mientras que el RMSE nos permite saber cuánto nos equivocamos en las unidades que medimos, en este caso en el número de espectadores al mes en los cines. En el Mean Error vemos que los primeros dos modelos suelen quedarse por encima del valor real, pues este es negativo, mientras que el del tercero es positivo, mostrando que suele quedarse por debajo del valor real.

AUTOMATIC SARIMAX-GARCH: MACHINE LEARNING APPROACH TOWARDS TIME SERIES PREDICTIONS

En este pequeño apartado vamos a contrastar los resutados que obtuvimos anteriormente con el ajuste automático, al que además le vamos a incluir una modelización de la varianza, para estimar los valores superiores e inferiores posibles.

SARIMAX_GARCH = function(d=datos,nendtrain=96, var="users", datecol="date"){
  
  preds = c()
  mean_val = c()
  results1 = c()
  results2 = c()
  
  for(i in seq(nendtrain, nrow(d)-12,12)) {
    
    train = data.frame(d[1:i, ])
    
    test = data.frame(d[(i+1):(i+12), ])
    
    train.ts = ts(train[,var], start=c(year(train[1, datecol]), month(train[1, datecol])), 
                   end=c(year(train[nrow(train), datecol]), month(train[nrow(train), datecol])),
                   frequency=12)
    
    test.ts = ts(test[, var], start=c(year(test[1, datecol]), month(test[1,datecol])),
                  end = c(year(test[nrow(test), datecol]), month(test[nrow(test),datecol])), 
                  frequency=12)
    
    xregs = as.matrix(calculoExplicativasCalendario(train$date,0)[,c("semanaSanta", "dt", "bisiesto", "festivo", "price")])

    
    ajusteAutomatico<-auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1, 
                             seasonal=TRUE,ic="aic",xreg=xregs,stepwise=TRUE)
    
    listaOutliers <- locate.outliers(ajusteAutomatico$residuals,
                                 pars = coefs2poly(ajusteAutomatico),
                                 types = c("AO", "LS", "TC"),cval=3)

   
    outliers <- outliers(c(as.character(listaOutliers$type)), c(listaOutliers$ind))
    outliersVariables <- outliers.effects(outliers, length(ajusteAutomatico$residuals))
    
    
    RegressorsMasOutliers <- cbind(xregs,outliersVariables)

    newxregs=as.matrix(calculoExplicativasCalendario(test$date, 0)[, c("semanaSanta", "dt", "bisiesto", "festivo", "price")])
    
    outliersNew <- outliersVariables[(i-11):i,]

    newxregs <- as.matrix(cbind(newxregs,outliersNew))
    
    #test = btc2$btc[i:i+5] #we set this to have 5 points so that the out.sample can work.
    m1.mean.model <- auto.arima(train.ts,max.d=1, max.D=1, max.p=1, max.q=1, max.P=1, max.Q=1, 
                             seasonal=TRUE,ic="aic",xreg=RegressorsMasOutliers,stepwise=TRUE)
    
    predicc = predict(m1.mean.model, n.ahead=12, newxreg= newxregs)
    
    preds = c(preds, predicc$pred)

    ar.comp <- arimaorder(m1.mean.model)[1]
    d_ <- arimaorder(m1.mean.model)[2]
    ma.comp <- arimaorder(m1.mean.model)[3]
    library(rugarch)
    model.garch = ugarchspec(mean.model=list(armaOrder=c(ar.comp,d_, ma.comp), arfima =F),
                             variance.model=list(garchOrder=c(1,1)),
                             distribution.model = "std")
    model.garch.fit = ugarchfit(data = train.ts, spec=model.garch, solver = 'lbfgs' )
    modelfor=ugarchforecast(model.garch.fit, data = NULL, n.ahead = 12)
    mean_val = c(mean_val, modelfor@forecast$seriesFor[,1])
    results1 <- c(results1, predicc$pred + modelfor@forecast$sigmaFor[,1])
    results2 <- c(results2, predicc$pred - modelfor@forecast$sigmaFor[,1])
  }
  
  real = data.frame(d)[(nendtrain+1):nrow(d), var]
  
  df = data.frame(pred=preds,meanval=mean_val,upper=results1, lower=results2, real=real)
  return(df)
}

Medidas de Error para el Modelo Automático

##                 ME     RMSE     MAE       MPE     MAPE
## Test set -729.9265 1650.171 1321.42 -11.12311 17.00692

Vemos que el ajuste automático deja unos resultados bastante pobres en comparación con el mejor modelo de los anteriores, que es el que hacemos con \(y^\lambda\). En aquél caso teníamos un 11.65 de MAPE, aquí tenemos un 17… Por lo tanto nos quedaríamos con el modelo 3 del apartado anterior, pues éste es el que mejor predice. Vamos a volver a dibujar el gráfico de las predicciones de una forma más completa, introduciendo también las predicciones del modelo automático.